import numpy as np
import matplotlib.pylab as plt
# 加载图像
im = plt.imread("BTD.jpg") # 加载当前文件夹中名为BTD.jpg的图片
print(im.shape) # 输出图像尺寸
# (4608, 2592, 3)即(y轴像素点数， x轴像素点数，图像通道数)
# 这里用的是RGB三通道图像，通道数为3

def plti(im, **kwargs):
    """
    画图的辅助函数
    """
    plt.imshow(im, interpolation="none", **kwargs)
    plt.axis('off') # 去掉坐标轴
    plt.show() # 弹窗显示图像

im = im[400:3800,:2000,:]  # 直接切片对图像进行裁剪
plti(im)

fig, axs = plt.subplots(nrows=1, ncols=3, figsize=(15,5)) 
# 将一张图分为1x3个子图，axs为各子图对象构成的列表。figsize为显示窗口的横纵比。

for c, ax in zip(range(3), axs): # 使用zip来同时循环3通道和3个子图对象
    tmp_im = np.zeros(im.shape) # 初始化一个和原图像大小相同的三维数组
    # 注意 tmp_im 仍然是三通道
    tmp_im[:,:,c] = im[:,:,c] # 只复制某一通道
    one_channel = im[:,:,c].flatten() # 索引该通道并展平至一维
    print("channel", c, " max = ", max(one_channel), "min = ", min(one_channel)) # 输出该通道最大最小的像素值
    ax.imshow(tmp_im) # 在子图上绘制
    ax.set_axis_off() # 去掉子图坐标轴
# 注意以上 tmp_im 采用的是切片复制
plt.show()

#输出：
#channel 0  max =  220 min =  11
#channel 1  max =  203 min =  10
#channel 2  max =  185 min =  0

# assert( np.log(np.e) == 1.0)
# np.log 即 ln() - 以e为底的对数函数

def do_normalise(im):
    return -np.log(1/((1 + im)/257) - 1)
# 预处理函数
# im中的像素值为 [0, 255] 闭区间， 则 (1+im) 为 [1, 256]
# 先做 (1+im)/257 操作将值归一化到 (0, 1) 开区间内
# 再使用 sigmoid函数 的反函数，效果见sigmod函数图像
# -np.log(1/((1 + 0)/257) - 1) = -5.5451774444795623
# -np.log(1/((1 + 255)/257) - 1) = 5.5451774444795623

def undo_normalise(im):
    return (1/(np.exp(-im) + 1) * 257 - 1).astype("uint8")
# 预处理函数的反函数
# 即先使用sigmod函数，再将值变换到(0, 257)区间再减1，通过astype保证值位于[0, 255]
# 关于 astype("uint8") ：
# np.array([-1]).astype("uint8") = array([255], dtype=uint8)
# np.array([256]).astype("uint8") = array([0], dtype=uint8)

def rotation_matrix(theta):
    """
    3D 旋转矩阵，围绕X轴旋转theta角
    """
    return np.c_[
        [1,0,0],
        [0,np.cos(theta),-np.sin(theta)],
        [0,np.sin(theta),np.cos(theta)]
    ]
# np.c_[ ] 将列表中的元素在第二维上拼接起来
# np.c_[[1,2],[3,4],[5,6]] =
# array([[1, 3, 5],
#        [2, 4, 6]])


im_normed = do_normalise(im)
im_rotated = np.einsum("ijk,lk->ijl", im_normed, rotation_matrix(np.pi))
# 利用爱因斯坦求和约定做矩阵乘法，实际上是将每个RGB像素点表示的三维空间点绕X轴（即红色通道轴）旋转180°。
im2 = undo_normalise(im_rotated)

plti(im2)

import matplotlib
print(matplotlib.matplotlib_fname())

from matplotlib.animation import FuncAnimation

fig, ax = plt.subplots(figsize=(5,8))

def update(i):
    im_normed = do_normalise(im)
    im_rotated = np.einsum("ijk,lk->ijl", im_normed, rotation_matrix(i * np.pi/10))
    im2 = undo_normalise(im_rotated)
    # 在更新函数里根据i来改变旋转的角度
    ax.imshow(im2)
    ax.set_title("Angle: {}*pi/10".format(i), fontsize=20)
    ax.set_axis_off()
    # 将旋转后的图绘出
# 以上其余注释见前文

#anim = FuncAnimation(fig, update, frames=np.arange(0, 20), interval=50)
# fig是图像句柄
# update是更新函数
# frames为帧数列表，将值依次提供给更新函数
# interval表示每帧间隔ms数

#anim.save('colour_rotation.gif', dpi=80, writer='imagemagick')
#plt.close()

def to_grayscale(im, weights = np.c_[0.2989, 0.5870, 0.1140]):
    """
    取原始图像的RGB值的加权平均来将图片转换为灰阶，权重矩阵为tile
    """
    # 默认的 weights = array([[ 0.2989,  0.587 ,  0.114 ]])
    tile = np.tile(weights, reps=(im.shape[0],im.shape[1],1))
    # assert( tile.shape == im.shape )
    return np.sum(tile * im, axis=2)
    # np.sum意味着沿某一轴求和，axis=2为第三维（0为第一维）
    # 整个乘法意味着由图像每个像素点的RGB 得到 (R*0.2989+ G*0.5870+ B*0.1140)灰阶值，图像的二维尺寸不变，而减为单通道。
img = to_grayscale(im)
plti(img, cmap='Greys') # 注意要以灰度形式画出

from scipy.ndimage.interpolation import zoom
im_small = zoom(im, (0.2,0.2,1))
# zoom 将图片每一维以相应系数缩小
# im.shape = (3400, 2000, 3)
# im_small.shape = (680, 400, 3)

from scipy.signal import convolve2d
# 引入二维卷积函数
def convolve_all_colours(im, window):
    """
    用窗口window卷积图像，依次对图像的每个通道卷积
    """
    ims = []
    # 用ims作为每个通道转换结果的暂存列表
    for d in range(3):
    # 对图像的三个通道循环处理
        im_conv_d = convolve2d(im[:,:,d], window, mode="same", boundary="symm")
        # mode决定输出尺寸，boundary决定边界条件，这里输出尺寸与原图相同，采用对称边界条件
        ims.append(im_conv_d)
        # 将单通道转换结果添加到列表

    im_conv = np.stack(ims, axis=2).astype("uint8")
    # 在第三维上堆叠ims列表中的每个元素，并通过astype保证值在0-255
    return im_conv

n=50
window = np.ones((n,n))
# 构建50x50的全1矩阵
window /= np.sum(window)
# 矩阵每个元素除以矩阵所有元素的和，使矩阵所有元素的和为1
plti(convolve_all_colours(im_small, window))

from scipy.ndimage import median_filter

def make_guassian_window(n, sigma=1):
    """
    使用高斯分布的权重创建一个n*n的方形窗口
    """
    nn = int((n-1)/2)
    a = np.asarray([[x**2 + y**2 for x in range(-nn,nn+1)] for y in range(-nn,nn+1)])
    # np.asarray可以将输入转化为np.array, 这里输入为一个列表推导式
    return np.exp(-a/(2*sigma**2))

def median_filter_all_colours(im_small, window_size):
    """
    对图像所有通道运用中值滤波
    """
    ims = []
    for d in range(3):
        im_conv_d = median_filter(im_small[:,:,d], size=(window_size,window_size))
        ims.append(im_conv_d)

    im_conv = np.stack(ims, axis=2).astype("uint8")
    
    return im_conv

window_sizes = [9,17,33,65]
fig, axs = plt.subplots(nrows=3, ncols=len(window_sizes), figsize=(15,15));

# 均值滤波 - 均匀窗口
for w, ax in zip(window_sizes, axs[0]):
    window = np.ones((w,w))
    window /= np.sum(window)
    ax.imshow(convolve_all_colours(im_small, window));
    ax.set_title("Mean Filter: window size: {}".format(w));
    ax.set_axis_off();
    
# 高斯滤波 - 高斯窗口
for w, ax in zip(window_sizes, axs[1]):
    window = make_guassian_window(w,sigma=w)
    window /= np.sum(window)
    ax.imshow(convolve_all_colours(im_small, window));
    ax.set_title("Guassian Filter: window size: {}".format(w));
    ax.set_axis_off();
    
# 中值滤波
for w, ax in zip(window_sizes, axs[2]):
    ax.imshow(median_filter_all_colours(im_small, w));
    ax.set_title("Median Filter: window size: {}".format(w));
    ax.set_axis_off();